% Compute the face space using some training face images
% training images are stored in ./trainfaces/*.pgm
%
% Author: Abhishek Dutta (adutta.np@gmail.com)
% June 04, 2008
%

%------------------------------------------------------
% STEP1: Read all 'M' training face images as Ti
%------------------------------------------------------

M = 10;     % no. of tranining faces
width = 92;
height = 112;
printf("\nReading training images ...");
T1 = pnmread("/root/octave/eigenfaces/trainfaces/1.pgm");
T2 = pnmread("/root/octave/eigenfaces/trainfaces/2.pgm");
T3 = pnmread("/root/octave/eigenfaces/trainfaces/3.pgm");
T4 = pnmread("/root/octave/eigenfaces/trainfaces/4.pgm");
T5 = pnmread("/root/octave/eigenfaces/trainfaces/5.pgm");
T6 = pnmread("/root/octave/eigenfaces/trainfaces/6.pgm");
T7 = pnmread("/root/octave/eigenfaces/trainfaces/7.pgm");
T8 = pnmread("/root/octave/eigenfaces/trainfaces/8.pgm");
T9 = pnmread("/root/octave/eigenfaces/trainfaces/9.pgm");
T10 = pnmread("/root/octave/eigenfaces/trainfaces/10.pgm");

%------------------------------------------------------
% STEP2: Compute average face vector Ψ
%------------------------------------------------------
%       Ψ = (1/M) Σ Ti      (i = 1, 2, ... , M)
% shi is the average face vector Ψ
printf("\nComputing average face vector ...");
shi = zeros(width,height);
% (T1, T2,...T10)/M not used as shi is a matrix that can only hold uint8 values
shi = (T1/M) + (T2/M) + (T3/M) + (T4/M) + (T5/M) + (T6/M) + (T7/M) + (T8/M) + (T9/M) + (T10/M);

%------------------------------------------------------
% STEP3: Compute the difference (Φi) of each face (Ti) from average face (Ψ)
%       Φi = Ti - Ψ
%------------------------------------------------------
% NOTE: @TODO phiX can be negative values ???
printf("\nComputing matrix A ...");
A = zeros(width*height,M);
A(:,1) = (T1 - shi);
A(:,2) = (T2 - shi);
A(:,3) = (T3 - shi);
A(:,4) = (T4 - shi);
A(:,5) = (T5 - shi);
A(:,6) = (T6 - shi);
A(:,7) = (T7 - shi);
A(:,8) = (T8 - shi);
A(:,9) = (T9 - shi);
A(:,10) = (T10 - shi);

%------------------------------------------------------
% STEP4: Computation of eigen vectors of covariance matrix (AA') is impractical.
%       Hence perform computation using (A'A).
%       Here,
%       A[N2xM] = [ Φ1 Φ2 .... Φm ], A[N2xM] and A'[MxN2]. Thus A'A[MxM]
%
%       STEP4.1:    Compute the Eigenvectors (Vi) of A'A matrix
%       STEP4.2:    the largest Eigenvectors (Ui) of AA' matrix is given by
%                   Ui = A * Vi
%       STEP4.3:    Normalize Ui such that || Ui || = 1 {ie: perform (1/||Ui||) * Ui }
%------------------------------------------------------

% DIM(m_cov_A) = MxM
printf("\nComputing A'A ...");
m_cov_A = A'*A;
% DIM(V_i) = M x 1
[V_i,mui_i] = eig(m_cov_A);
% DIM(U_i) = N^2 x 1
printf("\nComputing eigenvectors ...");
U_i = A * V_i;

% normalize U_i
% Verify correct computation of normal Ui using norm(U_i(:,1))
[row_U_i,col_U_i] = size(U_i);
printf("\nNormalizing eigenvectors ...");
for i = 1:col_U_i
    U_i(:,i) = U_i(:,i) / norm(U_i(:,i));
end

% represent faces using this basis
% reshape(A(:,j),height,width) is used to convert N^2x1 vector to a matrix of size NxN
for j = 1: col_U_i
    pnmwrite(strcat("eigfaces/eigen_face_",int2str(j),".pgm"),reshape(A(:,j),height,width), true)
end
printf("\nWritten eigenface images to eigfaces/*.pgm\n");    
